module space

  use globals
  
  implicit none
  
contains
!--------------------------------------------------------------------------

subroutine grid

  implicit none
  
  integer  :: ii, ecell
  real(8)  :: espace

  ecell = 7
  dr = rad/70.0d0   ! first cell size
  
  do ii = ecell, NI
	r(ii) = rad + dr*float(ii-1)*(1.0030d0**float(ii-1))
  enddo
  
  espace = r(ecell)
  
  do ii = 1, ecell-1
    r(ii) = rad + (espace-rad)/float(ecell-1)*float(ii-1)
  enddo
	  
end subroutine grid
!------------------------------------------------------------------------------
subroutine qdiff
  implicit none

  integer  :: ii
  
  Qr(1) = Qw_p
  Qr(NI) = 0.0d0
  
  do ii=2,NI-1
	 Qr(ii) = -cond*4.0d0*pi*r(ii)**2.0d0 * diff1(Tw,r,ii)
  enddo
    
end subroutine qdiff

!--------------------------------------------------------------------------
! First order finite difference for non-equal grids
function diff1(X,y,i)

  implicit none
  real(8) :: X(NI)
  real(8) :: y(NI)
  integer :: i
  real(8) :: diff1
  
	 diff1  = X(i+1)*(y(i)-y(i-1))/((y(i+1)-y(i-1))*(y(i+1)-y(i)))        &
                & + X(i)*(y(i+1)+y(i-1)-2*y(i))/((y(i+1)-y(i))*(y(i)-y(i-1))) &
		& + X(i-1)*(y(i+1)-y(i))/((y(i+1)-y(i-1))*(y(i-1)-y(i)))

end function diff1
!---------------------------------------------------------------------------
! Second order finite difference for non-equal grids
function diff2(X,y,i)

  implicit none
  real(8) :: X(NI)
  real(8) :: y(NI)
  integer :: i
  real(8) :: diff2
  
	 diff2  = 2*X(i+1)/((y(i+1)-y(i-1))*(y(i+1)-y(i))) &
                & + 2*X(i)/((y(i)-y(i+1))*(y(i)-y(i-1)))   &
                & + 2*X(i-1)/((y(i-1)-y(i))*(y(i-1)-y(i+1))) 

end function diff2
!------------------------------------------------------------------------------
! First order forward difference for equal grids
function fdiff1(X,y,i)

  implicit none
  real(8) :: X(NI)
  real(8) :: y(NI)
  integer :: i
  real(8) :: c0, c1, c2, c3, c4
  real(8) :: fdiff1
  
     c0 = -25.0d0/12.0d0
	 c1 = 4.0d0
	 c2 = -3.0d0
	 c3 = 4.0d0/3.0d0
	 c4 = -1.0d0/4.0d0
  
	 fdiff1 = (c0*X(i)+c1*X(i+1)+c2*X(i+2)+c3*X(i+3)+c4*X(i+4))/(y(i+1)-y(i))

end function fdiff1
!---------------------------------------------------------------------------
! Second order forward difference for equal grids
function fdiff2(X,y,i)

  implicit none
  real(8) :: X(NI)
  real(8) :: y(NI)
  integer :: i
  real(8) :: c0, c1, c2, c3, c4, c5
  real(8) :: fdiff2
  
	 c0 = 15.0d0/4.0d0
	 c1 = -77.0d0/6.0d0
	 c2 = 107.0d0/6.0d0
	 c3 = -13.0d0
	 c4 = 61.0d0/12.0d0
	 c5 = -5.0d0/6.0d0
  
	 fdiff2 = (c0*X(i)+c1*X(i+1)+c2*X(i+2)+c3*X(i+3)+c4*X(i+4)+c5*X(i+5))/((y(i+1)-y(i))**2.0)

end function fdiff2
!------------------------------------------------------------------------------
subroutine boundary(T,T_temp)
! wrong as is, not used.  TO BE UPDATED!
  implicit none
  real(8), dimension(NI) :: T, T_temp
  
  T(1) = Ta - Qw/(As*intCond)
  T(NI) = Tamb
  
  T_temp(1) = Ta - Qw/(As*intCond)
  T_temp(NI) = Tamb 
  
end subroutine boundary
!------------------------------------------------------------------------------

! numerical integration over the surface of a sphere using trapezoidal rule
! this assumes a full 2pi x pi amount of data is supplied
! Note: THIS IS FOR 3-PHOTON PROCESS ONLY !!!

function trap3(qdat,nx,ny)

  implicit none
  integer                      :: ii, jj, nx, ny
  real(8), dimension(nx,ny)    :: qdat
  real(8), dimension(nx)       :: temp
  real(8)                      :: trap3, w, rii, rnx, xphi, oo

  trap3 = 0.0d0
  temp = 0.0d0
  oo = 3.0d0
  rnx = nx

  do ii=1,nx ! polar
     rii=ii
     xphi = -90.0d0 + 180.0d0*(rii-1.0d0)/rnx
     w = abs(cos(xphi*pi/180.0d0))   ! weight to account for diminishing grid size towards poles
     do jj=1,ny-1 	             ! integral over each "ring", latitude / azimuth
	    temp(ii) = temp(ii) + (2.0d0*pi)/ny * w *(qdat(ii,jj)**oo+qdat(ii,jj+1)**oo)/2.0d0
     enddo
  enddo
  
  do ii=1,nx-1 ! polar
     ! integral over longitude
     trap3 = trap3 + (pi/nx) * (temp(ii)+temp(ii+1))/2.0d0
  enddo


end function trap3
!------------------------------------------------------------------------------

end module space

